Aggregate geometry in amyloid fibril nucleation 
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We present and study a minimal structure-based model for the self-assembly of peptides into 
ordered /3-sheet-rich fibrils. The peptides are represented by unit-length sticks on a cubic lattice 
and interact by hydrogen bonding and hydrophobicity forces. By Monte Carlo simulations with 
>10^ peptides, we show that fibril formation occurs with sigmoidal kinetics in the model. To 
determine the mechanism of fibril nucleation, we compute the joint distribution in length and width 
of the aggregates at equilibrium, using an efficient cluster move and flat-histogram techniques. This 
analysis, based on simulations with 256 peptides in which aggregates form and dissolve reversibly, 
shows that the main free-energy barriers that a nascent fibril has to overcome are associated with 
changes in width. 



Many proteins and peptides share the ability to self- 
assemble into amyloid fibrils, aggregates with a cross- 
P structure and remarkable mechanical properties, that 
are associated with a range of disorders as well as with 
functional roles The formation of amyloid fibrils, 

usually monitored by thiofiavin T (ThT) fiuorescence, is 
known to occur with reproducible sigmoidal kinetics [3], 
indicating a nucleation-dependent process. A powerful 
method for interpreting the experimental kinetic profiles 
is by means of rate equations This approach can re- 
veal some general properties of intermediate species par- 
ticipating in the growth process. It has proven useful 
for some related self-assembly phenomena as well, such 
as hemoglobin S aggregation f5^ and microtubule assem- 
bly [6 . Another method to elucidate the mechanisms of 
amyloid formation is by phase equilibria analysis [7]. 

By coarse-grained structure-based approaches [F, ad- 
ditional insights have been gained into the nucleation 
of amyloid fibrils [5l414|. Fibrillation pathways involve, 
however, a host of different aggregated species of widely 
varying size, and studying the competition among these 
species without restrictive assumptions represents a chal- 
lenge even in coarse-grained models. 

In this Letter, we introduce a minimal structure-based 
model that describes amyloid fibril formation in terms 
of physically inspired peptide-peptide interactions and 
yet allows for representative sampling of the model state 
space for relatively large systems. Using flat-histogram 
methods fl5l 116] and an efficient cluster move resembling 
the Swendsen-Wang algorithm for spin systems |17) . we 
determine equilibrium distributions in size and shape of 
the aggregated structures, in order to elucidate the free- 
energy landscape that a nascent fibril has to navigate. 

We consider N identical peptides, represented by unit- 
length sticks on a periodic cubic lattice with dimensions 

. We assume that the internal dynamics of a peptide 
are fast compared to the timescales for fibril formation, 
and therefore can be averaged out. 

Each peptide i is centered at a lattice site, Tj, and 




FIG. 1. Schematic illustration of the geometry and interac- 
tions of the model, (a) The orientation of a peptide is defined 
by the backbone vector (thick line) and the hydrogen-bond 
direction (thin line). The side-chain direction Si (dots) is 
given by = x p^. (b) Parallel /3-structure. (c) Antiparal- 
lel /3-structure. (d) Hydrophobic attraction. The /3-structure 
definitions are such that a /3-sheet has one hydrophobic and 
one polar side. 



two peptides cannot simultaneously occupy the same site. 
Associated with each peptide are two unit vectors 
and that can point in any of the six lattice directions 
(Fig. Ilk); hi represents the N-to-C backbone orientation, 
whereas ±Pj are the directions in which hydrogen bonds 
can form. The vectors and Pj are perpendicular, leav- 
ing a total of 24 possible orientations of a peptide. The 
vectors = b^ x Pj and — represent side-chain direc- 
tions. The +Si and — sides of a peptide are assumed 
to have different interaction properties and are referred 
to as hydrophobic and polar, respectively. 

The energy function describing the interactions be- 
tween the peptides is assumed pairwise additive, E = 
Si<j^iii where tij < 0. The pair potential is 
nonzero only if (i) peptides i and j are nearest neigh- 
bors on the lattice, and (ii) b^ and bj are perpendicular 
to r.ij — Tj — Ti and aligned either parallel or antiparallel 
to each other. When these conditions are met, we set 
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Eij = — 1 except in the three cases illustrated in Fig. [T] 
The first two cases correspond to parallel (Fig. [T]d) and 
antiparallel (Fig. [l];) /3-structure, respectively, and the 
third (Fig.[l]i) to hydrophobic side-chain attraction. The 
corresponding interaction energies are given by 



-(l + flp) 




parallel /3-structure 
antiparallel /3-structure (1) 
hydrophobic attraction 



The hydrophobic attraction is included because of evi- 
dence suggesting that a pairwise (steric zipper) /3-sheet 
organization is a common architecture for the core of 
amyloid fibrils [TB]. The b parameter must not be too 
large, in order for extended /3-sheets to form. Because 
the /3-sheets often are parallel in amyloid fibrils, we take 
Op > Oap, but the model can also be studied for flap > Op. 
In what follows, for simplicity, we stick to a single param- 
eter set, namely Op = 5, Cap = 3 and 6 = 1. With this 
choice, the parallel /3-strand organization dominates, but 
the suppression of antiparallel strand pairs is not pro- 
hibitively strong. 

We simulate the thermodynamics of this model using 
single-peptide as well as cluster moves. A cluster up- 
date makes it possible for aggregates to move without 
having to be first dissolved and then reassembled. To 
be able to also split and merge aggregates, we follow a 
stochastic Swendsen- Wang-type cluster construction pro- 
cedure [T7]. The construction is recursive and begins by 
picking a random first cluster member, i. Then, all pep- 
tides j interacting with peptide i {cij < 0) are identified, 
and added to the cluster with probability pij = 1 — e^'^'^ , 
where /3 — l/fceT is inverse temperature. This step is 
iterated until no cluster member has any unchecked in- 
teraction partner. Finally, the resulting cluster is subject 
to a trial rigid-body translation or rotation, drawn from 
a symmetric distribution, which is accepted whenever it 
does not cause any steric clashes. It can be verified that 
this algorithm fulfills detailed balance with respect to the 
canonical ensemble Pi, oc e~^-^" . 

To further enhance the sampling, we employ 
generalized-ensemble methods [TSl [16], along with 
reweighting techniques [19' . After estimating the den- 
sity of states, g{E), by the Wang-Landau method ^16j . 
we simulate the ensemble Pi, oc l/g{E^) |TS], where the 
distribution of E is flat. This approach was recently used 
for atomic-level aggregation simulations |20j and is use- 
ful for the present system as well, which displays phase 
coexistence at the fibrillation temperature, (see be- 
low). Our simulations sample a limited energy range, 
Emin < E < 0. The cutoff -Emin IS needed to prevent the 
formation of artificial cyclic aggregates, which otherwise 
may occur due to the periodic boundary conditions, but 
is sufficiently low to permit studies of temperatures in 
the fibrillar phase. 

The above cluster update can be adapted for the 
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FIG. 2. Overall thermodynamic properties of the A'^ = 256, 
L = 64 system at Tm ^ 0.6714. (a) Energy distribution. The 
sampled range is —1300 < -B < 0. Smoothing was applied to 
remove short-scale irregularities in g{E). (b) Mass fraction of 
aggregates with mass m, p{m), against m. Summed proba- 
bilities for three regions in m are indicated. Statistical errors 
on both p{E) and p{m) are comparable to the line width. 



generalized-ensemble simulations by adding an ac- 
cept/reject step, with acceptance probability Pacc(i^ 
v') = riivi\V,g{E^)Q-f^^'' lg{Ey.Y-f^^-'\. Here, /3 changes 
its meaning to become a tunable algorithm parameter. 
We did not fine-tune /3, but expect the optimal /3 to be 
in the neighborhood of /3„i = l/feBTm, as supported by 
preliminary runs. 

Using these methods, we studied the thermodynamics 
of the model for several different system sizes. Here, we 
focus on the results obtained for N = 256 and L = 64, 
corresponding to a peptide concentration of p ss 10"'^ per 
unit volume. This system size would have been very time- 
consuming to study with standard Monte Carlo methods. 

In our simulations, two distinct major phases occur: 
a high-energy phase dominated by small aggregates and 
a low-energy phase where large fibril-like aggregates are 
present. As displayed in Fig.|2^, at the midpoint temper- 
ature, Tin ss 0.6714, the energy distribution is bimodal, 
showing that the two phases coexist. Fig. [2]d shows the 
aggregate mass distribution, p{m), at T^^, which gives the 
probability for a peptide to be part of an aggregate with 
m peptides (m = 1 corresponds to free monomers). Like 
the energy distribution, p(m) is bimodal. The mass frac- 
tions of aggregates with to < 6, 6 < to < 62 and to > 62 
are 81.4%, 1.9% and 16.7%, respectively. Small aggre- 
gates are present in both phases and constitute a large 
fraction of the total mass at Tm (see also Supplemental 
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FIG. 3. Mass fraction of aggregates with length I and width 
w, p{l, w), at Tm for N = 256 and L = 64. 



Material, Fig. SI EB)- 

At first glance, the bimodality of p(m) may seem to in- 
dicate that fibril nucleation occurs when a critical mass 
is reached. However, this picture is geometrically incom- 
plete, because the species involved are neither strictly 
one-dimensional nor sharing one common shape, such as 
spherical. A simple but useful way to extend the anal- 
ysis is via the inertia tensor. As measures of the length 
and width of an aggregate, we define I = ^/12Xi + 1 and 
w = VI2A2 -I- 1, where Ai > A2 are eigenvalues of the 
inertia tensor. In our model, there is no interaction be- 
tween longitudinally aligned peptides to support growth 
in a third dimension. With these definitions, for a rect- 
angular aggregate, I and w are the numbers of peptide 
layers in the two directions. 

Fig. [sjshows the probability p{l, w) for a peptide to be 
part of an aggregate with length I and width w, at T^. 
Consistent with Fig. [2|3, p{l,w) is highest for aggregates 
with / small and w sa 1. Among larger aggregates, a clear 
preference can be seen for even over odd values of w, re- 
flecting a pairwise /3-sheet organization, although aggre- 
gates with six or more layers are severely constrained by 
finite-size effects. A second trend is that single-layer ag- 
gregates are shorter than two-layer ones, which in turn 
are shorter than those with four layers. We expect this 
trend to persist beyond the four-layer level if the system 
is sufficiently large. These overall features ofp{l,w) are 
likely to be quite robust, although the locations of the 
different maxima depend on both T and p. 

The shape of p{l, w) has implications for how fibrils 
nucleate and grow in the model. It suggests that the 
main free-energy barriers faced by a growing aggregate 
are associated with changes in width, and it must increase 
in width to be able to grow. 

Having examined the thermodynamics of the model, 
we now turn to the aggregation kinetics, studied using 
constant-temperature Monte Carlo dynamics. Because 
of evidence that amyloid growth occurs by monomer ad- 
dition 125], we here use single-peptide moves only. The 
simulations start from random initial conditions and the 
temperature is T = 0.66. 
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FIG. 4. Monte Carlo kinetics at T = 0.66. (a) Mass of the 
largest aggregate, mo(t), against time, t, in two representative 
runs with A'' — 256 and L = 64. Circles indicate where the 
width of the largest aggregate, uiq, switches from wo < 3.5 to 
Wo > 3.5. Time is given in sweeps, where one sweep consists 
of A'^ single-peptide updates, (b) Mass fraction of fibril-like 
aggregates, Xf{t), in 10 independent runs with = 131072, 
L = 512. Inset: log-log plot of the average over the 10 runs vs. 
t. The dashed line corresponds to a cubic growth, Xf{t) oc t^. 



We first stick to the system size N = 256 and L = 64, 
which is useful for examining the formation of individual 
fibrils. Fig. shows the mass of the largest aggregate, 
mo(t), against Monte Carlo time, t, in two representa- 
tive runs. Both runs exhibit an apparent waiting phase, 
before a large aggregate suddenly appears. Unlike aggre- 
gates occurring in the initial phase, this large aggregate 
is stable to dissolution. Near the jump in mass, a switch 
occurs in the width of the largest aggregate, wg- With 
a tiny fraction of exceptions, wq is below 3.5 before and 
above 3.5 after the switch point. Interestingly, as indi- 
cated in Fig. [4^, this switch in width occurs immediately 
before the sharp increase in mass. This suggests that the 
change in width is a critical event that renders the aggre- 
gate growth-competent. This finding matches perfectly 
with the shape of the p{l, w) distribution (Fig. |3|. 

The kinetics can also be studied for much larger sys- 
tems, which makes it possible to test in a direct manner 
whether or not the model captures the sigmoidal behavior 
observed experimentally. Fig. [4]d shows the time evolu- 
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tion of the mass fraction of fibril-like aggregates, X[{t), 
in 10 independent runs for N = 131072 and L = 512 
(same concentration as before). We define an aggregate 
as fibril-like if w > 3.5, w being the width, which ensures 
stability to dissolution. 

Comparison of these 10 runs (Fig. |4j3) shows that, for 
this system size, the kinetics are indeed sigmoidal and 
highly reproducible, as observed in bulk experiments ^Jj. 
At the end of the runs, the simulation box typically con- 
tains between 40-50 spontaneously formed fibrils (see 
Supplemental Material, Fig. S2 t21j), with an average 
length and width of I ^ 210 and w ~ 7, respectively. 
Inspection shows that the nucleation of new fibrils stops 
after roughly 10^ Monte Carlo sweeps. Existing fibrils 
continue to grow beyond that point, but eventually X{{t) 
levels off, due to monomer depletion. 

Our kinetic simulations, which do not include fragmen- 
tation events, may be compared to the classical Oosawa 
theory for homogeneous polymerization [23,. In particu- 
lar, this theory predicts the initial growth to be quadratic 
in time. Our data (Fig. |4]d) are well described by a cubic 
growth for small t, X{{t) cx t^. That the exponent ap- 
pears to be different than it is in the Oosowa theory is 
expected, because nucleation involves more than a single 
step in our model. 

In this Letter, we have presented a simple model for 
amyloid formation, where the nucleation of fibrils can 
be studied without any prior assumptions on the struc- 
ture of the aggregates involved. The formation of aggre- 
gated structures with a few /3-sheet layers has been ob- 
served in many previous simulations, also at the atomic 
level |24l [25] . Here, we have used systems much larger 
than in previous studies, to be able to examine the in- 
terplay between aggregate length and width in fibril nu- 
cleation. Our study shows that in this model the width 
of a growing aggregate plays a key role; to reach a given 
length, a minimum width is required, and to increase in 
width the aggregate has to overcome major free-energy 
barriers. Due to these barriers, the nucleation of a fib- 
ril occurs in distinct steps. The present study focused 
on the spontaneous aggregation of free peptides, but the 
model may also be useful for studying surface-catalyzed 
aggregation and the effects of a confining geometry. 

As in any model, simplicity is both a strength and a 
limitation. The final width of a growing fibril depends 
most likely on geometric factors left out in our model, 
such as twist. The question of what sets the final width 
is therefore beyond the scope of the present work. The 
assumption that internal degrees of freedom can be inte- 
grated out may be a good approximation for small fiex- 
ible peptides, but is clearly poorly justified for a folded 
protein that has to partially unfold before aggregation 
takes place. Our model further ignores any possible co- 
operativity of the interactions involved |261 [27] . In our 
model, aggregation is a highly cooperative process, al- 
though driven entirely by pairwise additive interactions. 



The simulations were performed at the LUNARC fa- 
cility, Lund University. 



[1 
[2 
[3 

[4; 

[5 
[6 
[7 
[8 
[9 
[10 

[11 
[12 
[13 
[14 
[15 
[16 

[17 
[18 

[19 
[20 
[21 

[22 

[23 
[24 

[25 

[26 

[27 



F. Chiti and C. M. Dobson, Annu. Rev. Biochem., 75, 
333 (2006). 

T. P. J. Knowles and M. J. Buehler, Nat. NanotechnoL, 
6, 469 (2011). 

E. Hellstrand, B. Boland, D. M. Walsh, and S. Linse, 
ACS Chem. Neurosci., 1, 13 (2010). 

T. P. J. Knowles, C. A. Waudby, G. L. Devlin, S. I. A. Co- 
hen, A. Aguzzi, M. Vendruscolo, E. M. Terentjev, M. E. 
Welland, and C. M. Dobson, Science, 326, 1533 (2009). 

F. A. Ferrone, J. Hofrichter, and W. Eaton, J. Mol. Biol., 
183, 611 (1985). 

H. Flyvbjerg, E. Jobs, and S. Leibler, Proc. Natl. Acad. 
Sci. USA, 93, 5975 (1996). 

J. D. Schmit, K. Ghosh, and K. Dill, Biophys. J., 100, 
450 (2011). 

C. Wu and J.-E. Shea, Curr. Opin. Struct. Biol., 21, 209 
(2011). 

R. Pellarin, E. Guarnera, and A. Caflisch, J. Mol. Biol., 
374, 917 (2007). 

S. Auer, C. M. Dobson, M. Vendruscolo, and A. Maritan, 

Phys. Rev. Lett., 101, 258101 (2008). 

M. S. Li, D. K. Klimov, J. E. Straub, and D. Thirumalai, 

J. Chem. Phys., 129, 175101 (2008). 

J. Zhang and M. Muthukumar, J. Chem. Phys., 130, 

035102 (2009). 

S. Auer and D. Kashchiev, Phys. Rev. Lett., 104, 168105 
(2010). 

N. S. Bieler, T. P. J. Knowles, D. Frenkel, and R. Vacha, 
PLoS Comput. Biol., 8, el002692 (2012). 

B. A. Berg and T. Neuhaus, Phys. Lett. B, 267, 249 
(1991). 

F. Wang and D. P. Landau, Phys. Rev. Lett., 86, 2050 
(2001). 

R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett., 58, 86 
(1987). 

M. R. Sawaya, S. Sambashivan, R. Nelson, M. L Ivanova, 
S. A. Sievers, M. L Apostol, M. J. Thompson, M. Bal- 
birnie, J. J. W. WiUzius, H. T. McFarlane, A. 0. Madsen, 

C. Riekel, and D. Eisenberg, Nature, 447, 453 (2007). 

A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett., 
63, 1195 (1989). 

S. M. Jonsson, S. Mohanty, and A. Irback, J. Chem. 
Phys., 135, 125102 (2011). 

See Supplemental Material at http:/ /link.aps.org/' sup- 
plemental/... for illustrations of p{E,m) and a typical 
final configuration from the kinetic runs. 
S. R. Collins, A. Douglass, R. D. Vale, and J. S. Weiss- 
man, PLoS Biol., 2, e321 (2004). 
F. Oosawa and M. Kasai, J. Mol. Biol., 4, 10 (1962). 

D. Li, S. Mohanty, A. Irback, and S. Huo, PLoS Comput. 
Biol., 4, el000238 (2008). 

M. Cheon, I. Chang, and C. K. HaU, Biophys. J., 101, 
2493 (2011). 

K. Tsemekhman, L. Goldschmidt, D. Eisenberg, and 
D. Baker, Protein Sci., 16, 761 (2007). 

B. Linse and S. Linse, Mol. BioSyst., 7, 2296 (2011). 



